Une méthode d'identification pour un système linéaire à retards 



A method of identification for a linear delays System. 

François Ollivier 
ALIEN, INRIA Futurs & LIX, UMR CNRS 7161 
École polytechnique 
91128 Palaiseau CEDEX, France 
f rancois . ollivier@lix . polytechnique . f r 

Saïd Moutaouakil et Brahim Sadik 
Département de Mathématiques 
Faculté des Sciences Semlalia 
B.P. 2390, Avenue Safi, Marrakech, Maroc 
{sadik, s .moutaouakil}@ucam. ac .ma 

2 février 2008 



Abstract 

We provide a class of methods for the identification of a linear System with delay of the shape x^ n ' (t) = 
^27=o a i x \t) + bu{t — h). They allow the simultaneous identification of the parameters and delay, the 
observation of x and its derivatives, knowing only generic input u and output x. They are robust to the 
noise. Used in continuous time, they may follow the évolution of slowly varying parameters and noise. 
The method may be generalized to Systems with two or more delays, e.g. x" (t) + ax(t — h\) = bu(t — h-z). 

Résumé On introduit une classe de méthodes d'identification et d'observation pour un système linéaire à 
retard de la forme x^ n \t) = Y27=o aiX (t) + bu(t — h). Celles-ci permettent l'identification simultanée du 
retard et des paramètres, l'observation de x et de ses dérivées et ne supposent que la connaissance de la 
sortie x et de l'entrée u, sans autre hypothèse que leur généricité. Elles se révèlent assez robustes au bruit. 
Utilisées en temps continu, elles permettent de suivre l'évolution de paramètres ou de retards lentement 
variables. La méthode peut être généralisée à des systèmes avec plusieurs retards, e.g. x"(t) +ax(t — h\) = 
bu(t — /12). 



Abridged English version 

Time delay Systems are known to take an important place in many nelds of application. We refer to [5] 
and the références therein for more détails on the subject. We propose a method for parameters and delay 
identification, related to the method introduced by Fliess and Sira-Ramirez [3] [5] (see also [T]). 

We write the System X)"=o a i x ^ (t) + bu(t) = 0, with a n = — 1. The output is y = x+w where w stands for 
the noise which is such that ^ w(t)cItJ / (Ti — T%) goes to when T% — T\ goes to infmity. The main idea 

is to use a family of functions /j such that fj(T\) = fj (T'a) = for k < n. Let I x j := /(T)x(r)dT. 
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Integrating by parts, we get I x {%)j. = m). So we may estimate the values of the coefficients ai and 

6 by solving the System YI^—q (^(—l) l ail^ ^(0^ + bl u j j = for j — I, . . . , m, by the mean squares method. 

We can estimate x and its derivatives in the same way, using functions <?j < j < n such that g^\T\) = 
for < k < n and gf ) (T 2 ) = for < k < j, with g [ f {T 2 ) ^ 0. 

In practice we have used fj(t) = (Ta — É) n+J e _A( ' T2_t \ the intégration being done between — oo and the 
current time. A good approximation of the intégrais is obtained by integrating the System J' x = x — A J Xj o 
and J' x j — jJx,j-i — ^Jx,j if j > 0, with initial conditions Jj(0) — 0: Jj(t) tends quickly to 4j ; _ n+1 for À 
great cnough. Numcrical simulations arc givcn with example 1. 

In section 3, we consider a delay System X)T=o a i x ^ % \t) + bu(t — h), with a n = — 1. We use the notation 
I xJtTu T a = Jt* /(( t - T i)/( T 2 ~ Tx))x{t)cLt, where / is such that /«(Q) = /«(l) = 0, for j < m with 
m > n. Let En^ T 3 )^c,ti dénote the équation 

n k 

^{n-TirtaJ^^+YsHTi-Tj-XjM^ = 0, 

i=0 1=0 

where the term 0(h min ^ k+1 ' m ') is neglected. Solving it by the mean squares method for generic x and u 
and a generic set S of couples (T\,T2), we get approximations âi of the coefficients and an approximation 
of delay equal to h = bi/b . We can then replace u by uy-(t) — u(t — h) in Er Tl ,T 2 ),x,u m order to get 
the improved approximation and iterate the process. Better précision could be achieved if / is such that 
f(t -ti)= J2 P k =o c k (h)f {k) (t), for example with f(t) = sin m (irt). 

In example 2, we use / = sin 2 on [0, n] and investigate the précision of the évaluation depending on 
the size of the noise. See Table and Table 1, which shows how our method may be improved by using a 
nonlinear least quare method, initialized with our data. 

In example 3 the method is generalized to the two delays System x"(t) + ax(t — h%) = buit — h-z). The 
standard déviation of our computed values are given in table 2 for différent noise level, stared values being 
improved by using a nonlinear least square method initialized with our results. 

In section 4, we adapt the method of section 1 to the delay situation, solving the System 

n 

E (t-^^yj") + + bJu hh , fj + 0(h 2 ) = 0, 

i=0 3 

where u?{t) = u(t — h) for \h — h\ -C 1. The delay évaluation is then h + bi/bo- For t > To, when bo does 
not vanish any more and the évaluation h(0) + 61/60 is assumed to be close enough, we takc h' = Xhbi/bo, 
so that h will converge to h. 

Example 3 shows a simulation with h = 0.5. Greater delays could be considered by changing the time 
scale. Example 4 considers the case of a slowly varying delay. Scilab simulation files are available at url [Ï3] . 

1 Introduction 

Les systèmes à retard sont l'objet de recherches actives, en raison de leur importance dans de nombreux 
champs d'applications. Nous renvoyons à l'article de synthèse [S] et aux références incluses pour plus de 
détails. L'identification simultanée du retard et des paramètres demeure un problème difficile, mais crucial 
pour la mise en œuvre de méthodes de contrôle appropriées. 

Nous proposons une méthode qui peut s'interpréter comme une variante de la méthode d'identification 
et d'observation due à Michel Fliess et Hebertt Sira-Ramirez [4j|5], et qui est susceptible de s'adapter à des 
systèmes à retard. L'article [T] fournit une application plus directe de cette approche. On pourra consulter, 
par exemple, [HJ 021 EU El E El H] pour un aperçu d'autres approches. 
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L'introduction d'un facteur d'oubli permet l'identification et l'observation en continu, ce qui est décrit 
dans la section 2 pour des systèmes sans retard et dans la section 4 pour des systèmes à retard. 

La section 3 décrit une méthode d'identification a posteriori. Une simulation permet d'évaluer sa précision 
en fonction du niveau de bruit. Nous avons utilisé à titre d'exemple des systèmes d'ordre 2, qui correspondent 
à une classe de modèles très répandus (voir [3])- Les calculs numériques ont été effectués sur Dell Latitute 
D400 avec Scilab 3.0, à partir de formules littérales calculées en Maple 9.5. 

2 Identification d'un système sans retard 

Pour alléger les formules, nous écrirons le système sous la forme : Yl7=o a i x (t) + bu(t) — 0, avec a n = —1. 
La sortie est y = x + w, où w désigne le bruit, que nous supposons tel que (^J^ w(T)d,T^ /(Ti — T 2 ) tend 
vers quand T 2 — Xi tend vers l'infini. L'idée de base est d'utiliser une famille de fonctions fj,j = l,...,m 
telles que /j fc) (T x ) = /j fc) (T 2 ) = pour k < n. 

On pose I x j :— J^ 2 f{T)x{r)dT. En intégrant par parties, on aura / (i) .(*,) = —I .(&+!), si k < n, et 
donc I x d) j. = On peut donc calculer la valeur des coefficients ai et b en résolvant le système 

n 

£((-l)W S)/ w)+ 6/^=0 

ï=o 3 

pour j = 1, . . . , m, si m > n. Pour m > n, ce système sera résolu par la méthode des moindres carrés. 

Si l'on se donne une famille gj, < j < n telle que g^ k \Ti) = pour < k < n et g^\T2) = pour 
< k < j, avec gf{T 2 ) £ 0, on obtient alors I sW = {-lfl w + Etof- 1 )'^ 1 ^"^)! 1 ' 1 ^). 
On estime x(T 2 ), . . . , a^" -1 - 1 (T 2 ) en résolvant le système 

E E ((-^"^^(raMra)^) + (-lîV^o +blu, 9j =0, 

i=0 \£=0 3 / 

pour j = 0, . . . , n — 1. Ces équations se déduisent immédiatement de a i^x( i ). gj + bl u>gj = 0. 

Par exemple, on peut prendre fj(t) = (T 2 - t) n+: > e" A ( T2 -*), et 5 3 = (T 2 - ^e _A ( T2_t ). Ces fonctions 
satisfont les hypothèses pour Ti = — oo. Si l'on pose J x = x — \J x ,o et J' x j = jJ x j-i — A J x j si j > 0, avec 
les conditions initiales Jj(0) = 0, Jj(T 2 ) tend vers I x , gj pour j = 1, . . . , n — 1 et vers I x ,fj- n+1 pour j > n, 
quand T 2 tend vers l'infini. La convergence est rapide pour A suffisament grand. On a g'j = Xgj — jgj-i et 

une relation similaire pour les fj. On a donc 4('),/ 3 = J2k=j-i Mi,j,klx,g k + J2l=i^i,n+j-i,klx,f k , où les 
matrices Mij^ satisfont des relations de récurrence simples, ce qui facilite les calculs. On obtient en un temps 
t générique et suffisament grand une bonne approximation des paramètres ai et b en résolvant le système 

n I n-l j \ 

E"' E M iJ:kIy,g k + E M i,n+j-l,k lyjk ' 1,1 - / = °i 
t=0 \fe=j-l k=l J 

que nous noterons C'A — D, où A = (ai, . . . , a n _i, 6)*. Comme il se peut que pour certaines valeurs le rang 
du système ne s'annule, il est préférable de ne pas le résoudre directement mais d'en chercher la solution 
par une méthode de gradient, en se donnant de nouvelles variables ôj et b, et en intégrant le système 
A' = —AC t (CA — D). Ce choix permet de conserver des valeurs précises lorsque C devient mal conditionnée 
et de réduire encore l'influence du bruit. Les valeurs de A et de A résultent d'un compromis entre la vitesse 
de convergence souhaitée et la précision recherchée. 
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Exemple 1 On a pris x solution de l'équation x" = a±x' + oqx + bu, où u = 60cos(l,23i + l,3sin(i) — 
0, 7cos(0, 5t)), avec les conditions initiales x(0) — 20 et x'(0) = 0,3. On a choisi un bruit gaussien avec un 
écart type de 5. La sortie est échantillonnée à 100 Hz. On a pris b = 2, ao = —0, 35, a\{t) = — 1, 2 si t < 30 
et ai = —1, 2 — 0, 02(t — 30) sinon. Les courbes de la figure^ résument les résultats obtenus pour À = 1 et 
A = 10~ 3 . On voit qu'après une phase de convergence, l'approximation obtenue est excellente. La précision 
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FiG. 1 - Estimations des coefficients de x et de x' , Estimations of the coefficients, x and x' . 



de l'évaluation des coefficients est naturellement moins bonne quand ai commence à varier. Toutefois, les 
évaluations de x et x' demeurent assez précises. 



3 Identification avec retards 

On considère maintenant un système à retard écrit sous la forme Yl7=o ^iX^ (t) + bu(t — h), avec a n = — 1. 
Soit / une fonction telle que /^(0) = /^(l) = 0, pour j < m avec m>n. On pose I x j,Tx,Tî = /(( T ~ 
T\) I ' (T2 — Ti))x(r)dT. En intégrant par parties, on déduit de l'équation du système 

]Tm - TJ-tatltjmnn + b ( £(T 2 - T x )-% Jit) <TuT A + (h^ k+1 ^) = 0. 
i=o \e=o / 

On désignera par E^ Tl T ^ x u cette équation, où le terme 0{h ral ^ k+1,m ^) a été négligé et chaque puissance h 1 
remplacée par bgjb afin de se ramener à une équation linéaire. Pour x et u génériques, on résoud par moindres 
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carrés le système E SiViU , s G S 1 en les inconnues ai et bt, où S est un ensemble générique d'intervalles d'un 
cardinal suffisant. On obtient des approximations a, et 60 des coefficients a» et b et une approximation du 
retard égale à bi/bo. La précision du résultat sera d'autant meilleure que le rapport h/(T 2 — Ti) est petit 
pour chaque couple de S. 

Une fois obtenue une approximation ho du retard, on améliore la précision du résultat en résolvant le 
système E SiViUho , s e S, où Uh (i) = u(t — ho), ce qui donnera la nouvelle approximation ho + b\/bo- On itère 
le processus. 

Pour de meilleurs résultats, on peut également choisir une fonction / telle que f(t + h) — X^fc=o a k{h)f^ k \ 
par exemple /(*) = t m (l - t) m ou sin m (7rf). 

Exemple 2 Nous avons fait des simulations numériques avec f(t) = sin 2 (t) sur l'intervalle [0, 7r] pour un 
système d'ordre 2. L'équation E/ TltT2 \ yu s'écrit alors (7\ — T2)~ 1 a\L x j^T 1 ,T 2 +aoIx,f,T 1 ,T 2 + boI u j,T 1 ,T 2 
-hI u ,cos(2t),T 1 ,T 2 +b2l u ,sin(2t),T 1 ,T 2 + 0(/i 2 ) = (Ti - T 2 )~ 2 I x j» >Ti ,t 2 et l'on obtient une approximation de h 
égale à ho + acos(&2/&o) à partir d'un système E s ^ Uho , s e S, où h est une estimation préalable du retard. 

Les valeurs des coefficients sont les mêmes que celles de l'exemple 1. On a choisi la commande u = 
60 cos(l, 23i + 0, 33 * sin(t) — 0, 47 cos(0, 5i)) et un retard h = 4. Le bruit est un bruit gaussien et la fréquence 
d'échantillonnage est de 500 Hz. On a pris pour S les couples (10k, 10k + 15), 1 < k < 9. On part de ho = 
et l'on itère les calculs une cinquantaine de fois (en un temps de l'ordre de quelques secondes sur un Dell 
Latitude D400), ce qui est suffisant pour atteindre un point fixe à la précision numérique près. 

Le tableau ci-dessous, qui donne la moyenne (M.) et l'écart type (E.) de séries de cent estimations des 
coefficients et du retard réalisées pour des bruits gaussiens d'écarts types croissant, montre que l'on obtient 
des estimations crédibles, même pour un niveau de bruit notable. 



E. bruit 


M. ai 


É. ai 


M. a 


E. a 


M. 6 


E. b 


M. h 


E. h 


1 


-1,20 


0,020 


-0,35 


0,007 


2,00 


0,04 


4,00 


0,015 


2 


-1,20 


0,028 


-0,35 


0,010 


2,00 


0,05 


4,00 


0,023 


5 


-1,17 


0,084 


-0.34 


0,029 


1,94 


0,17 


3,97 


0,065 


10 


-1,12 


0,15 


-0,32 


0,055 


1,84 


0,33 


3,91 


0,132 



Tableau 0. Moyenne et écart type des valeurs calculées 
Table 0. Mean and standard déviations (E.) of computed values 



Nous avons comparé nos résultats avec une méthode souvent utilisée en pratique qui consiste à évaluer les 
coefficients par une méthode de moindres carrés non linéaire, pour différentes valeurs du retard régulièrement 
espacées, et à choisir celle pour lesquelles on approche le mieux la sortie. Cette méthode est inapplicable pour 
l'exemple 2 avec la méthode de moindre carré non linéaire disponible dans scilab, car celle-ci ne converge vers 
les paramètres que si elle est initialisée avec des valeurs voisines. Toutefois, les moindres carrés non-linéaires 
permettent d'améliorer la précision de l'évaluation des paramètres et du retard, une fois initialisés avec les 
résultats fournis par notre méthode. Le tableau ci-dessous fournit les valeurs approchées des écarts types 
obtenus par notre méthode et après utilisation des moindres carrés (*), évaluées sur 100 essais, la durée 
moyenne d'un calcul étant de 10, 5 s. 



E. bruit 


E. ai 


É. ai* 


É. a 


É. a * 


E. b 


É. 6* 


É. ie(0) 


É. :r(0)* 


É. rr'(0) 


É. x'(0)" 


É. h 


E. h' 


5 


0,1 


0,002 


0,03 


0,0005 


0,2 


0,002 


1 


0,4 


6 


0,5 


0,07 


0,001 


10 


0,15 


0,005 


0,06 


0,001 


0,35 


0,005 


1,7 


0,6 


10 


0,9 


0,1 


0,0015 



Tableau 1. Écart type des valeurs calculées 
Table 1. Standard déviations of computed values, stared values improved with raean-square 



Ce type de calcul suppose, naturellement l'identifiabilité du système. Si la commande est nulle, la sortie 
constante etc. celui-ci ne sera pas identifiable et il ne le sera que localement pour un comportement périodique. 
Notons que, nos calculs reposant uniquement sur des résolutions au sens des moindres carrés, il n'y a jamais 
d'erreurs numériques du type « division par 0. » Pour des retard grands, les itérations peuvent diverger et 
produire une erreur du type invalid index si les intervales d'intégration sortent de la période d'observation 
du système. 

Exemple 3 Cette méthode se généralise au cas de plusieurs retards. Considérons le système x" (t) + ax(t — 
h\) = bu(t — /12). En utilisant f(t) = 1 — cos(t) sur [0,27r], on a aL u j.T 1 .T 2 +b sin(/i2/(T2 — Ti))/a;,sin,Ti,T2 
+0(hj) -bI xJ , Tl ,T 2 -asm(hi/(T 2 - Ti))/ XiS i niTl ,T 2 +0{h\) = (T 2 - T^Ixj- ,t x ,t 2 - En choisissant un 
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nombre suffisant d'intervalles [T\,T 2 ], on peut évaluer par moindres carrés linéaires des approximations de 
a, asin(/ii/(T 2 — T\)), b et bsin(hi/(T 2 — Ti)). On en déduit des valeurs de h\ et h 2 que l'on peut réinjecter 
dans les équations pour améliorer la précision. 

Les simulations ont été faites avec x = 3sin(t/2) + 2cos(i/3) ; a = 2, 7 b = 1, 5, h\ = 2, h 2 = 4. 
L'echantillonage est effectué à 500 Hz avec des bruits gaussiens d'écart-types différents. On a pris T 2 — T\ = 
10 et 20 intervalles avec T\ variant par pas de 2 à partir de Xi = 15. Le tableau suivant donne la moyenne et 
l'écart-type des valeurs calculées (sur 100 essais, chaque calcul durant environ 1,5 s. avec 10 itérations) pour 
différents niveaux de bruit. Les valeurs de h\ et h 2 ne sont naturellement définies que modulo la période 12ir 
de la sortie. La précision du calcul se dégrade si les intervalles n'ont pas une longueur T 2 — T\ suffisament 
grande devant h\ et h 2 , mais aussi pour un signal périodique si T 2 — T\ n'est pas assez petit devant la période. 



E. du bruit 


M. a 


E. a 


M. b 


E. b 


M. ht 


E. h t 


M. h- 2 


E. h 2 


0,025 


2,69 


0,012 


1,50 


0,006 


1,99 


0,005 


3,99 


0,0051 


0,05 


2,69 


0,018 


1,49 


0,010 


1,99 


0,010 


3,99 


0,011 


0,1 


2,70 


0,043 


1,50 


0,025 


1,99 


0,01 


3,99 


0,017 


0,2 


2,68 


0,095 


1,48 


0,054 


1,99 


0,039 


4,00 


0,043 



Tableau 2. Moyenne et écart type des valeurs calculées 
Table 2. Mean and standard déviations (E.) of computed values 



4 Identification en temps réel avec retard lentement variable 

On peut également adapter aux systèmes à retard les méthodes d'identification et d'observation en temps 
réel introduites dans la section 1. On utilise les équations 




avec bi = bh e et = u(t — h). Ceci suppose que \h — h\ <C 1. On peut s'y ramener quel que soit le retard 

par un changement d'échelle de temps. L'évaluation du retard est alors h + bi/b . On résoud le système 
par intégration au moyen d'une méthode de gradient comme décrit section 2. Lorsque t > T , on suppose 
que les valeurs des estimations sont devenues assez stables (il faut au moins que &o n e s'annule plus). On 
pose alors h' = si t < T et h' — Xhbi/b sinon. Pour plus de précision, on calcule I u „ f j en posant 

^jrfi + A^i/MKW-A/,,/,) 

Exemple 4 On utilise les mêmes commande, coefficients, type de bruit et fréquence que pour l'exemple 1. On 
prend T = 14, A = 1, A = 4.10~ 3 et \ h = 0, 4. 



Exemple 5 On utilise les mêmes commande, coefficients, type de bruit et fréquence que pour les exemples 1 
et 3. On prend T = 40, À = 1, A = 4.10~ 4 et Xh = 0,4. Le retard est initialisé à 3,6 sa valeur à T , afin 
d'accélérer la convergence. 

5 Conclusion 

Les méthodes que nous avons introduites permettent un large choix dans leur mise en œuvre. Les quelques 
essais réalisés montrent qu'elles sont applicables dans des situations réalistes, en particulier avec un niveau de 
bruit conséquent. Elles peuvent se généraliser aisément, au moins en théorie à des situations plus complexes : 
retards multiples ou n'intervenant pas uniquemement au niveau de la commande, comme l'illustre l'un de 
nos exemples. 



G 



1.19r 
0.88 




FiG. 2 - Estimations des coefficients et du retard, Estimations of coefficients and delay. 




FiG. 3 - Estimations du retard de x et de x', Estimations of the delay, of x and of x' . 



Leur utilisation suppose la connaissance a priori d'une borne supérieure pour le retard, ce qui est essentiel 
pour un choix approprié des différents paramètres A, A, À/j et Tq. En revanche, leur capacité à s'adapter à 
des coefficients ou à un retard lentement variable nous semble un atout intéressant. Les fichiers utilisés pour 
les simulations sont disponibles via l'url [T3] . 
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